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Theory of space-charge waves on gradient-profile relativistic electron beam: 

an analysis in propagating eigenmodes 

Gianluca GeloniQ Evgeni Saldin, Evgeni Schneidmiller, and Mikhail YurkmQ 

Deutsches Elektronen-Synchrotron, 
Notkestrasse 85, 22607 Hamburg, 
Germany 
(Dated: February 2, 2008) 

We developed an exact analytical treatment for space-charge waves within a relativistic electron 
beam in terms of (self-reproducing) propagating eigenmodes. This result is of obvious theoretical 
relevance as it constitutes one of the few exact solution for the evolution of charged particles under 
the action of self-interactions. It also has clear numerical applications in particle accelerator physics 
where it can be used as the first self-consistent benchmark for space-charge simulation programs. 
Today our work is of practical relevance in FEL technology in relation with all those schemes where 
an optically modulated electron beam is needed and with the study of longitudinal space-charge 
instabilities in magnetic bunch compressors. 



PACS numbers: 52.35.-g, 41.75.-i 
I. INTRODUCTION 

The evolution problem for a collection of charges under 
the action of their own fields when certain initial condi- 
tions are given is, in general, a formidable one. Numerical 
methods are often the only way to obtain an approximate 
solution, while there are only a few cases in which finding 
an exact treatment is possible. 

In this paper we report a fully self-consistent solution 
to one of these problems, namely the evolution of a rel- 
ativistic electron beam under the action of its own fields 
in the (longitudinal) direction of motion. The problem 
of longitudinal space-charge oscillations has been, so far, 
solved only from an electrodynamical viewpoint [1| , or us- 
ing limited one-dimensional models [2| . On the contrary, 
in our derivation the beam, which is assumed infinitely 
long in the longitudinal direction, is accounted for any 
given radial dependence of the particle distribution func- 
tion. 

An initial condition is set so that the beam is consid- 
ered initially modulated in energy and density at a given 
wavelength. When the amplitude of the modulation is 
small enough the evolution equation can be linearized. 
An exact solution can be found in terms of an expansion 
in (self-reproducing) propagating eigenmodes. 

Our findings are, in first instance, of theoretical impor- 
tance since they constitute one of the few exact solutions 
known up to date to the problem of particles evolving 
under the action of their own fields. 

Next to theoretical and academical interest of our 
study we want to emphasize here its current relevance 



'Electronic address: gianluca.aldo.geloni@desy.de Also at Depart- 
ment of Applied Physics, Technische Universitcit Eindhoven, The 
Netherlands 

tAlso at Joint Institute for Nuclear Research, Dubna, Moscow Re- 
gion, Russia 



to applied physics and technology. For example, particle 
accelerator physics in general and FEL physics in partic- 
ular make large use of simulation codes (see for instance 
0,0,13) m order to obtain the influence of space-charge 
fields on the beam behavior. Yet, these codes are bench- 
marked against exact solutions of the electrodynamical 
problem alone (i.e. solutions of Poisson equation) and 
only recently [2| partial attempts have been made to 
benchmark them against some analytical model account- 
ing for the system evolution. However, such attempts 
are based on one-dimensional theory which can only give 
some incomplete result. On the contrary, we claim that 
our findings can be used as a standard benchmark for 
any space-charge code from now on. 

To give some up-to-date example of practical applica- 
tions (besides, again, theoretical and numerical impor- 
tance) we wish to underline that our results are of rel- 
evance to an entire class of problems arising in state- 
of-the-art FEL technology. In fact, several applications 
rely on feeding (optically) modulated electron beams into 
an FEL. For instance, optical seeding is a common tech- 
nique for harmonic generation || . Moreover pump-probe 
schemes have been proposed that couple the laser pulse 
out of an X-ray FEL with an optical laser. To solve 
the problem of synchronization between the two lasers 
a method has been proposed, which makes use of a sin- 
gle optically modulated electron beam in order to gener- 
ate both pulses. This proposal relies on the passage of 
this electron beam through an X-ray FEL and an opti- 
cally tuned FEL pfl : given the parameters of the system, 
plasma oscillations turn out to be a relevant effect to 
be accounted for. It is also worth to mention, as an- 
other example, the relevance of plasma oscillation theory 
in the understanding of practical issues like longitudinal 
space-charge instabilities in high-brightness linear accel- 
erators. High-frequency components of the bunch current 
spectrum (at wavelengths much shorter than the bunch 
length) can induce, through self-interaction, energy mod- 
ulation within the electron beam which is then converted 
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into density modulation when the beam passes through 
a magnetic compression chicane, thus leading to beam 
microbunching and break up. When longitudinal space- 
charge is the self-interaction driving the instability (see 
0), one is interested to know with the best possible ac- 
curacy how plasma oscillations modulate the bunch, in 
energy and density, before the compression chicane. 

Our calculations can be applied to these cases directly 
or in support to macroparticle simulations thus providing 
outcomes of immediate practical importance. 

Throughout this paper we will make use of cgs units. 
Our work is organized as follows. Next to this Introduc- 
tion, in Section [HI we pose the problem in the form of 
an integro-differential equation which makes use of nat- 
urally normalized quantities. In Section IIIII we present 
our main theoretical results. In the following Section llVl 
we describe some applications and important exemplifi- 
cations of the obtained results, including the role of the 
initial condition. Finally, in Section we come to con- 
clusions. 



in energy and density is customary. Once the modula- 
tion wavelength X m has been fixed, it is natural to define 
the phase tp — w m {z/v z {£q) — t), where v z (£q) ~ c is 
the longitudinal electron velocity at the nominal beam 
kinetic energy £q — (7 — 1)toc 2 , u) m = 2irv z /A m , t is the 
time and z the longitudinal abscissa. Upon this, and af- 
ter what has been said about the choice of longitudinal 
dynamics, it is appropriate to operate in energy-phase 
variables (P, tp), P being the deviation from the nominal 
energy. 

In this spirit the total derivative of the phase tp is given 
by: 

dtp _ dtp dtp dt _ u m uj m 
dz dz dt dz v z (£q) v z (£) 

Now, if we assume that the particle energy is not signifi- 
cantly different from the nominal energy we can expand 
v(£) in £ around £q. Keeping up to second order terms 
in £ and using the definition of P we find 



II. 



THEORY 



A. Position of the problem 

We are interested in developing a theory to describe 
longitudinal plasma waves in a relativistic electron beam. 
In order to do so we have to find a way to translate in 
mathematical terms the idea of dealing with the longi- 
tudinal dynamics only, while keeping intact the general 
three-dimensional description of the system (electromag- 
netic fields and particle distribution). 

Immediately related with the development of a theory, 
is its range of applicability. The idea of selecting longi- 
tudinal dynamics translates from a physical viewpoint in 
the assumption that the longitudinal space-charge fields, 
alone, describe the system evolution. This corresponds 
to a situation with physical parameters tuned in such 
a way that, on the time-scale of a longitudinal plasma 
oscillation, transverse dynamics does not play a role. 

One may, of course, devise different methods to select 
the longitudinal dynamics alone: actually he will come 
up with some ideal approximation of what real focus- 
ing systems in modern linear accelerators do, provided 
that the /^-functions are large enough with respect to the 
plasma wavelength. 

Whatever the practical or ideal mean chosen to select 
longitudinal dynamics such a choice translates, from a 
mathematical viewpoint, in the choice of dynamical vari- 
ables along the direction of motion only, while transverse 
coordinates enter purely as parameters in the description 
of the fields and of the particle distribution. 

Our beam is initially modulated at some wavelength 
A m , in density and energy. This is no restriction because, 
as said in Section [IJ the modulation amplitude is consid- 
ered small enough so that we can linearize the evolution 
equation. Then, a Fourier analysis of any perturbation 



dtp 
dz 



(2) 



where we took advantage of the fact that 
(dv z /d£) \s = s ~ 0/(7^0), where j z = (1 - 
v z (£o) 2 1 'c 2 ) -1 / 2 . Note that, here, we distinguish 
from the very beginning between 7 and j z (or v and 
v z ). In fact our theory can be applied to the case of an 
electron beam in vacuum as well as to the case of a beam 
under the action of external electromagnetic fields, for 
example in an undulator: in the first situation 7 = j g , 
strictly, while in the latter they obviously have different 
values. 

The full derivative of P is simply given by 



dP 

dz 



= -eE z 



(3) 



where E z (z,tp) is the space-charge field in the z direc- 
tion. Eq. @ and Eq. @ are the equation of motion 
for our system and they can be interpreted as Hamil- 
ton canonical equations corresponding to the Hamilto- 
nian H(tp, P, z): 



H 



dtPE z 



2c 7 2 £0 



(4) 



In this sense, Eq. (0J, alone, defines our theory. The 
bunch density distribution will be then represented by 
the density / = f(tp, P, z; r±), where the semicolon sep- 
arates dynamical variables from parameters and it will 
be subjected to the (Vlasov) evolution equation 



df 



dH df 
~dPlhp 



dH df 
lhp~dP 



(5) 
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with appropriate initial condition at z = 0. Here we are 
interested in a beam initially modulated in energy and 
density; moreover, as said in Section HJ the modulation 
must be small enough to ensure that linearization of Eq. 
(0 is possible. Thus we will take f(ip,P,z;r±)^ =0 = 
fo(P;r±) + fi(ip,P,z;r±)\ z=0 . f is a so called unper- 
turbed solution of the evolution equation Eq. JSJ, and 
it does not depend on z, while /i is known as the per- 
turbation; it is understood that, in order to be in the 
linear regime, fi <C f for any value of dynamical vari- 
ables or parameters. In the following we assume that the 
dynamical variable P and the parameter r± are separa- 
ble in /o so that we may write fo(P] r±) = n (r±)F(P), 
where the local energy spread function F(P) is consid- 
ered normalized to unity. The initial modulation can be 
written as a sum of density and energy modulation terms: 
/i(V>, P, z; r ± ) lz=0 = hdW, P) r±) + fi e (i>, P; r±). 

On the one hand fid is responsible for a pure density 
modulation and can be written as 



Ji 



-ec 



dP ai d F + a lt 



-ej / dz' 
Jo 



E, 



dF 
dP 
.dF 



e ii, e »^3T 



1^0 



dP^-e^> 
dP 



(z'-z) 



(10) 



The next step is to present the equation for the electric 
field E z which, coupled with Eq. (|10|l . will describe the 
system evolution in a self-consistent way. 

We start with the inhomogeneous Maxwell equation 
for the z-component of the electric field 



„ 2 „ I d 2 E z A dp e 4tt dj z 

V 2 E Z = 47T— -| — , 

c 2 dt 2 dz c 2 dt 



(11) 



where p e is the electron charge density. Remembering the 
definition of complex quantities E z and ji and accounting 
for the fact that p e ~ j z /v z one can rewrite Eq. as 



P, z- n) = a u (r ± )F(P) cos(^) , 



(6) 



where we set to zero an unessential, initial modulation 
phase. 

On the other f\ e is responsible for a pure energy mod- 
ulation and can be assumed to be 



dF 

/i e (V>, P, z; r±) = a le (r±)—cos(ip + ip ) , (7) 



V 2 



d 2 E^ 1 d 2 E z e^ 



dz 2 



dt 2 



Andj^ And^e 



dz 



dt 



(12) 



where is the Laplacian operator over transverse coor- 
dinates. Explicit calculations of partial derivatives with 
respect to t and z give 



where ipo is an initial (relative) phase between den- 
sity and energy modulation. Finally it is convenient 
to define complex quantities fid — a±dF, and f\ e = 
a le (dF/dP)e 1 ^ so that f 1]z=0 = {fi d + fu)e^+CC. In 
the linear regime, then, one can write f\(ip,P,z;r±_) = 
/i(P,^;rj.)e^ + CC. Further definition of E z = 
E(z;r ± ) in such a way that E z = E z e^ + E* z e~ 1 ^ al- 
lows one to write down the Vlasov equation, Eq. JSJ, 
linearized in f\: 



— i — - 

dz cj 2 £ Q 



(8) 



Eq. (JHJ) is far from being the final form of the evolution 
equation, since we still have to couple it with Maxwell 
equations, which constitute the electrodynamical part of 
the problem. However an integral representation of f\ 
can be given at this stage: 



h = h 



c->* £ o +en — 
dP 



dz ' E z e 



-(z'-z) 



(9) 



Let us now introduce the longitudinal current density 
jz{z;r ± ) = -j (r x ) + ~iie 1 ^ + j*e-^, where j {r±) ~ 



ecno(r±) and j\ ~ — ec J_ 
grated in P thus giving 



dPfi . Eq. 10 can be inte- 



d 2 E z e^' 
dz 2 



d 2 E z 



2i 



uJ m dE z 
v z dz 



'■-E 7 



d 2 E z e^ 
dt 2 



= -JLeJ^ 



djie 



dz 



dji 

dz 



. UJ 7n ~. 

v z 



djie 



Zljj 



dt 



-lUJ m Jie 



, (13) 



(14) 



(15) 



(16) 



Substitution back into Eq. (|12|l yields (v z ~ c): 



V 2 ± E Z 



d 2 E z 



2i 



dE z io 2 m E z 



lie 2 



dz 2 v z dz 

47r dji Airiuj m ~ 

v z dz 



Iz 1 - 



It is reasonable to assume that the envelope of fields and 
currents vary slowly enough over the z coordinate, in 
order to neglect first and second derivatives with respect 
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to z in Eq. i|17|) . Mathematically, this corresponds to the 
requirements 



guided by Eq. 121(1 to introduce the transverse size pa- 
rameter 



dE z 



dz 



d 2 E z 



and 



dz 2 



dji 
dz 



E. 



E 7 



h 

™m |~ I 

7i 



(18) 



(19) 



(20) 



where we introduced the wave number k m = 2tt/ X m . We 
note that 2j z X m is, roughly, the field formation length: 
by imposing conditions l|18l) , (|19fl and (|20ll we are requir- 
ing that the characteristic lengths of variation for current, 
field and its derivative are much longer than the field for- 
mation length (actually condition i|19|) over the variation 
of the field derivative already included in condition <|18|) . 
since 7^ 3> 1): this simply means that we can neglect 
retardation effects or, in other words, that the fields are 
known at a certain time, when the charge distribution are 
known at the same time. This assumption is not a re- 
striction and it is verified in all cases of practical interest. 
Then Eq. l(T7|) is simplified to 



q = k m r /j z 



(23) 



Eq. I|21|) can now be written in its final dimensionless 
form: 



V\E Z - q 2 E z 



i<l 2 ji 



(24) 



where is the Laplacian operator with respect to nor- 
malized transverse coordinates. 

By analyzing Eq. I|1(J|I and using the normalized quan- 
tities defined above we recover a dimensionless variable 
z = Apz, where Ap is given by 



A P =[4//(^ 77z 2 )] 1/2 



(25) 



I A = mc 3 /e being the Alfven current. 

Using Eq. (|25|l and looking now at the exponential 
factors in Eq. (|1(J|I it is straightforward to introduce the 
dimensionless energy deviation P = P/(p£o), p being 
defined by 



(26) 



V\E Z 



lie 2 



47HW m - 
2 2 Jl 

lie 2 



(21) 



which forms, together with Eq. (|10f> . a self-consistent 
description for our system. 

Similarity techniques can be now used in order to ob- 
tain a dimensionless version of Eq. (|10|l and Eq. (|21|l . 
First note that the dependence of jo on the transverse 
coordinates can be expressed, in all generality, as 



jo = IoS {r ± /r ) 



S (r±/r )dr ± 



(22) 



where Io is the beam current, tq the transverse pro- 
file parameter (i.e. the typical transverse size of the 
beam), So the transverse profile function of the beam and 
the integral in dr± is calculated over all the transverse 
plane. Furthermore it is understood that the normaliza- 
tion of Sq is chosen is such a way that So(0) = 1. It 
is then customary to introduce the current density pa- 
rameter J = Io [J S(r±/r )dr±] so that we can de- 
fine quite naturally the dimensionless current densities 
jo = jo/ Jo = S {r±/r ) and ji = ji/J - It follows 
from Eq. (|21|l that the electric field should be normal- 
ized to Eq = AnJo/LUrm which suggests the definition 
E z = E z /Eq. For the normalization of the transverse 
coordinates we use f = r±/rQ so that one is naturally 



From the definition of P it follows immediately that the 
factor p£o is a natural measure for energy deviations. 
The rms energy spread ((A£) 2 ) can be measured by the 
dimensionless parameter 



((AS) 2 ) 



n 2 f 2 
P c 



(27) 



The local energy spread distribution F(P) was de- 
fined as normalized to unity, so that it is customary to 
introduce F(P) as the distribution function in the re- 
duced momentum P also normalized to unity. For ex- 
ample, when the energy spread is a gaussian we have 
F(P) = (27rA|)^ 1 /2 e ^P 2 /(2A 2 r ). When A| < 1 our 

beam can be considered cold meaning that F(P) ~ S(P). 
However note that, in order to specify quantitatively the 
range of validity of the cold beam assumption with re- 
spect to the values of Aip, one should first solve the more 
generic evolution problem for a non-cold beam and then 
study the limit for small values of K\. 

Now Eq. I|10|) can be expressed in the final, dimension- 
less form: 



.71 



dP d ld F + a lc 



dF 

'dP 



-iPz 



6 



-So / dz' 
Jo 



E- I <ir—r r ■' 

dP 



(28) 



In the case of a cold beam F(P) 
transforms to 



S(P) and Eq. 



where did = —ecaid/Jo and d\ e = — ece"^°ai e /( JqpEq). 
One can combine Eq. I|28[l and Eq. I|24|) in order to 
obtain a single integrodifferential equation for E z or, al- 
ternatively, an integral equation for j\ . 

As regards the description of the evolution in terms of 
Eg, direct substitution of Eq. (J2EJ m Eq. fZQ yields 
immediately 



E x 



-i {d\d + id\ e z) + Sq dz' (z! — z) E 



(33) 



Finally, after double differentiation with respect to z 
we get back the well-known pendulum equation for one- 
dimensional systems: 



V] E - <fE, = Iq 2 I dP !d u F + a le ^ j e- iPi 



-iq 2 So / dz' 
Jo 



E 7 



, dP 



(29) 



On the other hand, the description of our system in terms 
of ji can be obtained first by solving Eq. (|24|l and then 
substituting the solution in Eq. (|10|l . For the solution of 
Eq. H24fl we can use the following result (see for example 
0): 



E, 



n 

' 2vr 



dr^nKo (q | 



~(s) 



(30) 



where Kq indicates the modified Bessel function of the 
second kind. Then, substitution in Eq. (|2"5|) yields 



d 2 E z 
dz 2 



SqE z 



. 



(34) 



In this approximation, and in the particular case Sq — 
1, we recover the one-dimensional plasma wave number 



A p , which agrees with the normalization 



A, 



It 



should be noted that such a normalization is natural in 
the limit q — > oo but it progressively loses its physical 
meaning as q becomes smaller and smaller: of course, 
using our dimensionless equations for q 1 will still 
yield correct results, because the equations are correct, 
but the normalization fits no more the physical feature 
of the system, in that case. 

Eq. I|32l) . which describes the system evolution in the 
limit q — > oo, corresponds to the the one-dimensional 
case. We can use Eq. (13211 and impose £ = thus getting 
the electromagnetic field at the beginning of the evolu- 
tion, i.e. the solution of the electromagnetic problem in 
the limit q — > oo: 



Ji = 



dP ( d ld F + ai e ^C ) e" 
. I dP 



iPz 



d pdF e iP iS '- 2) 
dP 



(31) 



Note that the description of the system in terms of fields 
or currents is completely equivalent. Using one or the 
other is only a matter of convenience; it will turn out 
in the following Sections that the description in terms 
of the fields is particularly suitable for analytical ma- 
nipulations, while the description in terms of currents is 
advisable in case of a numerical approach. 

It is interesting to explore the asymptote of Eq. I|29|) 
for q — > oo. In this case one obtains the following simpli- 
fied equation for the field evolution: 



E, 



dP [ di d F + die^C ) e 
\ dP 



-iPz 



iSr, / dz' 



E z 



00 dP^JPiz'-Z) 

-oo dP 



(32) 



Ez = -iji U=0 ■ (35) 
This can be easily written in dimensional form as 



E 7 



4? 



h 

c 



(36) 



where j\ ~ Ii/(7rr§). 

We can check Eq. (|36|l with already known results in 
scientific literature. In fact, the field generated by an 
electron beam modulated in density in free space can be 
easil y c alculated in the system rest frame (see e.g. 
and In the limit of a pancake beam (i.e. for large 

transverse dimension) the following impedance per unit 
length of drift has been found: 



Ai 



(37) 



The first factor on the right hand side of Eq. 1)36(1 , which 
is the impedance per unit drift according to our calcula- 
tions, is exactly the result in Eq. 1|37|) which proves that 
our starting equation Eq. l(2"9"|l can be used to solve cor- 
rectly the electromagnetic problem in the case q — > oo, 
as it must be. 



7 



Finally, before proceeding, it is worth to estimate the 
value of our parameters for some practical example and 
to see how our assumptions compare with an interest- 
ing case. Nowadays photo-injected LINACs, to be used 
as linear colliders or FEL injectors constitute cutting 
edge technology as regards electron particle accelera- 
tors. Currents of about 3 fcA with energy such that 
7 ~ 10 3 can be achieved together with quite a small 
energy spread {(AS) 2 ) 1 / 2 ~ 100 fceV, radial dimensions 
of order r ~ 100 /im and a normalized emittance e n ~ 1 
mm mrad. In these conditions, a typical modulation 
wavelength of about 1 /im will lead to a value of q ~ 1, 
which is well outside the region of applicability of the 
one-dimensional theory so that a generalized theory like 
the ours must be used. It is interesting to note that, 
by means of the definitions of q and A p , we can write 
p = [I/(jlA)] 1 ^ 2 q~ 1 ' this suggests that a natural mea- 
sure for the current is indeed jIa which is . in practice, 
the Alfven-Lawson current (see 01 > E3)- I n our 
numerical example q ~ 1 and the value of p is actually 
determined by the ratio I /(jIa) ~ 2 • 10~ 4 . As a result 
p ~ 10~ 2 : this means that our simplified Maxwell equa- 
tion, Eq. (|21jl . is valid up to an accuracy of 10~ 2 . More- 
over, in this case, A 2 - ~ 10~ 4 so that the cold beam case 
turns out to be of great practical interest. Finally e„ ~ 1 
mm mrad. This corresponds, for 7 ~ 10 3 , to a betatron 
function (3f ~ 77q/ c n ~ 10 m. Our choice of consider- 
ing the longitudinal motion alone is satisfied when the 
period of a plasma oscillation is much shorter than the 
period of a betatron oscillation, that is when A p /3/ 3> 1. 
It should be noted, though, that this estimation cannot 
have a rigorous mathematical background before a more 
comprehensive theory, including the effects of transverse 
dynamics, is developed: only then our present theory can 
be reduced to an asymptote of a more general situation, 
and conditions for its applicability can be derived in a 
rigorous way. Keeping this fact in mind, in our example 
A p ~ 3 • 10" 1 m~ 1 which means A p f3f ~ 3 signifying that 
this particular case is at the boundary of the region of 
applicability of our theory. 



III. 



MAIN RESULT 



Given Eq. (|29[l with appropriate initial conditions for 
it is possible to find an analytical solution to the 
evolution problem. The method is similar to the one 
used for the solution of the self-consistent problem in 
FEL theory (see and relies on the introduction of 
the Laplace transform of E z , namely: 



following ordinary differential equation: 

V\-q 2 {l-iDS a )\E = iq 2 (b Q a ld + Da le ) . (8!)) 
where 



Do = / dP- 



F 



and 



00 p + iP 



D=l dP d -^t 

-00 p + iP 



(40) 



(41) 



with the boundary conditions E — ► for |fj_| — ► 00 
and dE/dr± — > for |fj_| — > 00. We can rewrite Eq. 
Pty as 



CE = f 



having introduced 



C = Vl + g(r±,p) 



f(*x,p) = iq 2 [D a ld + Den 



and 



g(v ± ,p) = -q 2 (l-iDS ) 



(42) 



(43) 



(44) 



(45) 



Note that only Eq. I|29|l can benefit from the use of the 
Laplace transform but not the integral equation Eq. 13111 . 

Eq. 112(1 is a nonhomogeneous, linear, second-order 
differential equation. We are interested in solving Eq. 
(|4*2")l for any given p such that Re(p) > 0. Solution is 
found if we can find the inverse of the operator C, namely 
a Green function G obeying the given boundary condi- 
tions; in this case we simply have 



E= / dr' ± G(r ± ,r' ± )f(r' ± ) 



(46) 



A. Generic approach 



E(p,t±) = / dze~v ? 'E z , (38) 
Jo 

with Re(p) > 0. The advantage of the Laplace transform 
technique is that the evolution equation is transformed 
from the integrodifferential equation Eq. I|29|) into the 



Depending on the choice of g, i.e. on the choice of So, 
F and p, the differential operator C can change its char- 
acter completely making £ more or less difficult to deal 
with. For example, the case of a self-adjoint operator is 
obviously a simple situation, since its eigenvalues are real 
and its eigenfunctions form a complete and orthonormal 
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set for the space of squared-integrable functions L 2 (de- 
fined over the entire transverse plane through fj_) with 
respect to the internal product: 



<h\g>= J dr ± gh* . (47) 

Then, such a orthonormal set can be used to provide, 
quite naturally, an expansion for G. 

However, in the most general situation, C is not self- 
adjoint: to see this, it is sufficient to note that g is not 
real. As a result, the eigenvalues of C are not real, its 
eigenfunctions are not orthogonal with respect to the in- 
ternal product in Eq. (|47l) . we do not know wether the 
spectrum of C is discrete, completeness is not granted and 
we cannot prove the existence of a set of eigenfunctions 
either. 

To the best of our knowledge there is no theoretical 
mean to really deal with our problem in full generality. 
When not self-adjoint operators are encountered in differ- 
ent branches of Physics (see, for example, 0] and 0) 
mathematical rigorousness is somehow relaxed assuming, 
rather than proving, certain properties of the operator. 
We will do the same here assuming, to begin, the exis- 
tence of eigenfunction sets; then, as for example has been 
remarked in 0] and 0] one can consider, together with 
the spectrum of C defined by the eigenvalue problem: 

V&i = A,*, (48) 

also the spectrum of its adjoint, defined by 

= [A,-]*** . (49) 

It can be shown by using the bi-orthogonality theorem 
[H that 



< **|^ >= / dPjtfj-tfi = S jt . (50) 

In other words the sequence admits {^jjj as a 

bi-orthonormal sequence. Then one has to assume com- 
pleteness and discreteness of the spectrum so that the 
following expansion is correct: 

G = J2\G*j><**\, (51) 

3 

We ascribe to alternative theoretical approaches and nu- 
merical techniques the assessment of the validity region 
of this assumption, which should be ultimately formu- 
lated in terms of a restriction on the possible choices of 
So and F. In other words we give here a general method 
for solving our problem which is valid only under the ful- 
fillment of certain assumptions, but we make clear that 
it is not possible, to the best of our knowledge, to strictly 



formulate the applicability region of this method in terms 
of properties of Sq and F as it would be desirable. 

With this in mind one can use the fact that G = C^ 1 
and write 

C(f± , fl) = E 2#i&21). m 

3 3 

Finally, substituting Eq. in Eq. lj4^|) one gets 

E = E ^Xf 1 J ■ (53) 

To find E z we use the inverse Laplace transformation 
that is the Fourier-Mellin integral: 

E z (z,v ± ) = — dpE(p,i x )f?* , (54) 

^ 7r * J a—ioo 

where the integration path in the complex p-plane is par- 
allel to the imaginary axis and the real constant a is pos- 
itive and larger than all the real parts of the singularities 
of E. 

The application of the Fourier-Mellin formula comes 
with another, separate mathematical problem related 
with the ability of performing the integral in Eq. I|54|) . 
One method to calculate the integral is to use numerical 
techniques and integrate directly over the path defined, 
on the complex p-plane, by Re(jp) — a. 

Yet, there is some room for application of analytical 
techniques left. In fact, under the hypothesis that E is 
also defined, except for isolated singularities, as an ana- 
lytical function on the left half complex p-plane and on 
the imaginary axis and under the hypotesis that E — > 
uniformly faster than l/|p| fc for a chosen k > and for 
Arg(p) within [ir/2, Sit/2] one could use Jordan lemma 
and close the integration contour of Eq. (|54|) by a semi- 
circle at infinity on the left half complex p-plane. An 
obvious (and well-known) problem is that E is defined 
only for Re(p) > according to Eq. (J3SJ. Yet, if the 
border points at Re(p) = are regular points of E (ex- 
cept for isolated singularities) then one can consider the 
(unique) analytical continuation of E along the border, 
from the original domain of analyticity (i.e. the points 
p with Re(p) > except for isolated singularities) to 
the entire complex plane (again, isolated singularity ex- 
cluded). Then one can still apply Jordan lemma on the 
analytic continuation of E (provided that it obeys the 
other assumption), because the final result is the inte- 
gral in Eq. (|54H which is uniquely defined by the original 
function E for Re(p) > 0. 

The problem is trivially solved for the case of a cold 
beam because F = S(P) so that Dq = 1/p and D = 
—i/p 2 . Then Eq. I|53|) defines indeed an analytic function 
in all points of the complex plane with the exception 
of p = and the points such that Aj(p) — 0. All the 



9 



hypothesis of Jordan lemma are verified and the method 
can be applied without any problem. 

The situation is completely different in the case of a 
generic energy spread function F. In fact by inspec- 
tion of Eq. i|40[) and Eq. (|41() one is immediately con- 
fronted with the fact that the integrands in Dq and D 
are, usually, singular at all the points of the imaginary 
axis Re{p) — since the integration in P is taken from 
—oo to +oo. As a result the points Re{jp) = are not reg- 
ular points of E and E cannot be analytically continued 
through the border Re(p) = 0. 

This problem is the same encountered in the treatment 
of Landau damping (see Of course one may follow 

the solution proposed by Landau and present particular 
definitions of D and D at Re(p) = that are 



f°° - F 

D Q = (P) / dP r + nF(ip) Re{p) = (55) 

J -oo p + iP 

and 

D = (P) r dP dF/dP + nF'(ip) Re(p) = , (56) 
J-oo p + iP 

so that Do and D are now regular for Re(p) = and can 
be (uniquely) continued at Re(p) < by 



D = I dP + 2-Kptip) Re(p) < (57) 

> p + iP 



and 



D = 



dP 



dP/dP 
P + iP 



+ 2TiP\ip) Re(p) < . (58) 



In this way the definition of E could be extended, ex- 
cept for isolated singularities, to an analytical function 
on the entire complex p-plane and, if F(P) behaves rela- 
tively well, Jordan lemma can be applied without further 
problems. 

Yet, we think that the application of Landau's pre- 
scription, i.e. the definition in Eq. (|55() and Eq. I|56l) 
should be taken with extreme caution. As it is reviewed 
by Klimontovich [TEI ] and references therein, Landau's 
method is equivalent to the introduction of additional as- 
sumptions on the system, namely the adiabatic switching 
on of the space-charge field at t = — oo. 

Another method equivalent to Landau's consists, as 
has been remarked long ago by Lifshitz |T^ |. in the in- 
troduction of a small dissipative term into the linearized 
Vlasov equation which ceases to be non-dissipative from 
the very beginning. The Vlasov equation is then solved 
by Fourier technique and the limit for a vanishing dissi- 
pating term is taken in the final result, which leads, in 
the end, to Landau's result. Yet, the limit for a vanishing 
dissipation must be taken in the final formulas in order to 
recover Landau's coefficient and not before. This means 



that Landau's method consists in the introduction of ad- 
ditional assumptions regarding the system under study 
or equivalently, in changing the very nature of the equa- 
tions describing our system (from non-dissipative to dis- 
sipative): therefore in practical calculations, we prefer to 
deal only with the case of a cold beam where the orig- 
inal non-dissipative nature of the system is maintained 
without problems, leaving the other case to future study. 
Such a viewpoint constitutes a restriction but as we have 
seen in the previous Section the cold beam case is prac- 
tically quite an important issue: in this Section, we will 
present our results in full generality without fixing F but 
keeping in mind, however, all the warnings discussed be- 
fore. 

In any case and again, in all generality, we can say 
that wether or not the conditions of Jordan lemma are 
satisfied depends on the distribution P(P). In the case 
Jordan lemma is applicable one can find a closed, analytic 
expression for P z : 



E z (z,r^ = J2t>A?±y 



dp 



/•OO 

/ dfl$„(f' ± )f(f' ± ,Aj) , (59) 
Jo 



where $>j(r±) = ^j(r±,p = Xj) and Xj are solutions of 
the equations: 



A i (p) = 



(60) 



or, which is the same, solution of Eq. (|48|l as Aj — 
0: from this viewpoint the functions Q n j constitute the 
kernel of the operator L and Xj are the values of p such 
that C admits a non-empty kernel. 

It is interesting to note that { } j is a subset of {^j}j, P 
naturally suited to expand any function of physical in- 
terest (the field E z ). In this sense, one may say that the 
fields are subjected to constraints given by Maxwell equa- 
tions, which are codified through the operator £; these 
constraints are implicitly used during the anti-Laplace 
transform process, thus selecting only those {^j}j :P of 
physical interest. In this spirit, although p is mathemat- 
ically allowed to span over all the complex plane with 
Re(p) > 0, only the particular values for which p = Xj 
have physical meaning in the final result. 

An explicit expression for (dhj(p) / dp) p= \ j can be 
found. Using Eq. I|5l)|l and Eq. (|48|l we can write 



(61) 



Aj= JdrxVjWj ■ 
Then, differentiating Eq. I|61|) we have 



dp 



op J 



(62) 
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Our final result is therefore written as follows: 



which may be rewritten using Eq. (|66|l as: 



where the coupling factor Uj is given by 



(63) 



J*i(g) * 



(64) 



Eq. (|63[) describes the evolution of the system under the 
action of self-fields in a generic way, for any bunch trans- 
verse shape So j f° r anv choice of local energy spread P 
and any initial condition (under the assumptions men- 
tioned before). Our solution is indeed an analysis in 
(self-reproducing) propagating eigenmodes of the electric 
field. 

We have seen that, due to the fact that C is in general 
not self-adjoint, the modes ^ are not orthogonal in the 
sense of Eq. (|4*T|) nor, as a consequence, $j are. More- 
over, even if ^> j were orthogonal, $j are chosen among 
the tyj at different values of p so that orthogonality of 
<j>j with respect to Eq. I|47|) is also not granted, ft is 
possible, however, to formulate appropriate initial condi- 
tions to obtain a single propagating mode as a solution 
of our self-consistent problem. This demonstrates that 
single modes have physical meaning besides being math- 
ematical tools for function decompositions. 

Suppose, for example, that we wish to excite a single 
mode at fixed values of j. On the one hand Eq. (|63[) is 
simplified to 



OE z 



A, 



au 



(70) 



Finally, if a le ^ 0, comparison of Eq. J7UJ) and Eq. (|(77|) 

gives: 



au 



(71) 



Since for plasma oscillations Xj has to be imaginary, Eq. 
(|71|l fixes the phase ipo = ^ (with I integer number). 
The actual shape of did (and au) is obtained, modulus 
a multiplicative constant, by substitution of Eq. I|65|l. 
calculated at z = 0, in Eq. 16 7|) and it is fixed by the 
following condition: 



aid = 0(^0 



(72) 



We will now present some remarkable example of how 
to apply Eq. i|63[l and explore, in particular cases, the 
applicability region of our method. 

We will start our exploration discussing the situation 
of an axis-symmetric beam which is still quite a generic 
one. Given the symmetry of the problem we will make 
use, from now on, of a cylindrical (normalized) coordinate 
system (f, 0, z), with obvious meaning of symbols. 

Since ji = ji(f, z, 0), E z = E z (f,z,4>) and / = 
f(r,p,(j)), it is convenient to decompose them in az- 
imuthal harmonics according to 



E z = Uj $j ef 



(65) 



and, differentiating with respect to z, one also obtains 



dE z 
dz 



XjE z 



(66) 



On the other hand, the evolution equation, Eq. (|29f) at 
z = reads: 



OE z = iq aid , 

I z— 



where we introduced 



O = VI - q 2 



(67) 



(68) 



The same Eq. 11291 differentiated with respect to z and 
evaluated at z — gives 



ji{z;f,4>) 



E z (z;r,cb)= £ E^\z;r)e 



in<i> 



and 



f(r,p,<f>)= J2 f {n) (r,p)e- m ' p ■ 

n— — oo 

Moreover, in cylindrical coordinates we have 



(73) 



(74) 



(75) 



ld_ 

f df 



Of 



+ 



1 d 2 



(76) 



These definitions allow us to write our equations and 
results for the n-th azimuthal harmonic of the electric 
field. In this situation the operators C and O can be 
written as 
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where, now, So = So(f) in the definition of g and 



O 



so that 



d 2 1 d 
dr 2 f dr 



(78) 



Within this special situation we will now treat in de- 
tail the case of a stepped or parabolic transverse profile. 
Further on we will see how the solution for the stepped 
transverse profile can be used to obtain a semi-analytic 
solution for any transverse profile. 



(n) 



(79) 



Having specialized our results to the axis-symmetric case 
it is worth to spend some words on the nature of the oper- 
ator C in one particular case. As we have said at the be- 
ginning of this section, the case of a self-adjoint operator 
is a particularly blessed one. It is interesting to note that 
in the axis-symmetric case, when So, F and p are such 
that g is real, we deal with a singular Sturm-Liouville 
problem as it is shown immediately by multiplying side 
by side by f Eq. J75J). In fact, in this case C can be 
written in the usual form presented in Sturm-Liouville 
theory: 



dr 



dr 



- — + rg(f,p) 



(80) 



In this case given the internal (axisisymmetric) product 
inL 2 : 



<h\ g>~- 



dffgh* 



(81) 



self-adjointness condition for C is satisfied in the interval 
[0, oo) by all the functions in the linear space § which 
we define as the space of integrable-square functions not 
singular with their first derivatives at r = and such 
that, for any /, g chosen in § the following condition is 
satisfied: 



lim f[9*(r)f'(f)-f(f)g'*(f)} = 



(82) 



Note that § is a subset of L 2 . 

This is of course a very particular situation, interesting 
to discuss but unfortunately not very useful in practise, 
since in order to solve our problem we still have to as- 
sume that Eq. (|51|l is correct for a generic p. When 
this assumption is made, in the axis-symmetric case our 
results Eq. (J^Sl and Eq. take the form 



B. Stepped profile 

Consider the case of a step function So = 1 for f < 1 
and So = for f > 1. In this case g(f,p) inside the 
operator C is simply given by 



-r ii 

„2 



f < 1 
f > 1 



(85) 



Let us restrict to the assumption of a cold beam with 
F(P) = S(P). Then D = i/p 2 . 

First we look for the solutions of Eq. I|48|l with the 
boundary condition that ip n j and their first derivatives 
vanish at infinity. The search for the eigenfunctions can 
be broken down into an internal (f < 1) and an exter- 
nal (r > 1) problem, with the conditions of continuity for 
ip n j and its derivative across the boundary, since the final 
result, the electric field, is endowed with these properties 
too. We recognize immediately that the internal and the 
external problems are, respectively, the complex Bessel 
and modified Bessel equations with appropriate bound- 
ary conditions. 

Keeping in mind the physical nature of our problem, 
we impose that ip n j must be regular functions of f over 
[0, oo). Then, without loss of generality, we can exclude 
Bessel functions Y n and /„ from entering our expression 
for 1p n j. 

Since for the field calculations we are interested in find- 
ing the eigenfunctions 4> n j = ip n j(r,p — v ' 

impose = 0, thus obtaining 



Xj""') we can 



C 2 K n (qf) 



where a, 



Aj(p) = 0, to be still determined at this stage. Im- 
posing continuity of ip n j and its derivative at f = 1 one 
finds 



f < 1 
f > 1 

An) 



(86) 



i(r,p)\ and are roots of 



,j nj ^nj 



(r)e> 



where 



u nj (f) = 



f™d?r'$ nj fW(r',\f) 



(83) 



(84) 



C*2 — C\ 



Jn{0ij) 

K n (q) ' 



(87) 



which leaves the choice of an unessential multiplicative 
constant, and 



ajJ' n {a 3 )K n {q) - qK' n (q)J n (a 3 ) = 



(88) 
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Eq. (|88|l can be rewritten with the help of recurrence 
relations for Bessel functions in the following form: 

a.jJ n+ x(aj)K n (q) - qK n+1 (q)J n (aj) = . (89) 

Eq. (|89|l is our eigenvalue equation, defining the val- 
ues of aj or, equivalently, of A^. Since q is real and 
positive one must have that a.jJ n +i(otj)/ J n (ctj) is real 
and positive. Then, it can be shown that aj must 
be real. As a result \j are imaginary and such that 
-1 < Jm(A^ n) ) < 1. For any given eigenvalue , also 

— Xj is solution of Eq. lj%9l corresponding to fast and 
slow plasma waves: from now on we will consider, for sim- 
plicity of notation, only the branch Im(X^) > 0. Note 

(n) 

that from a physical viewpoint, the condition that Xj is 
imaginary means that we are in the absence of damped 
or amplified oscillations. On the other hand, the fact 
that Im{xf Y ) < 1 means that plasma oscillations have a 
minimum wavelength given by 2tt/A p . 

It is interesting to plot the behavior of Im{X^), pa- 
rameterized for several values of n and j as a function 
of q. Fig. Fig. and Fig. 03 show the behav- 
ior, as a function of q, of the first five eigenvalues for 
the first three azimuthal harmonics. It should be noted 
that Im(\j ) increases with q and therefore with ro, 

but A p scales as r^ 1 ; as a result, the period of the self- 
reproducing solution identified by fixed values of n and j, 

that is 27r/(A p /m(Aj"' 1 )), will increase as ro is increased. 
As it can be seen by inspection all the imaginary parts 
of the eigenvalues converge to 1 as q — ► oo; this can 
also be derived directly from Eq. (JS5J. As q — ► oo 
we have K n (q) ~ K n+ i(q), so that Eq. lj%9l gives sim- 
ply J n {otj) — 0; this is possible only when aj — f n ,j, 
where v n j is the j-th root of J n . Then, in this limit, 

(xf) 2 = -q 2 /(q 2 + ^j) and (a^)* -> -1 since 

q — > oo. Note that convergence to unity tends to get 
slower as n and j increase. On the other hand, when q 
becomes smaller and smaller the plasma wavelength as- 
sociated with each mode starts to differ significantly from 
A p and, as noted before, we should use an effective A p in 
place of A p in our dimcnsionless quantities in order for 
these to retain their physical insight. 

The asymptotic behavior of 7m(A) for 5 « 1 can be 
derived directly from Eq. (|89|l too. We consider first the 
case n = 0. In the limit g<lwe have qKi(q) j Ko(q) ~ 
— (In q — In 2 + 7^) , where je is the Euler gamma con- 
stant. We remember that xJi(x)/Jq ~ x 2 /2 for x 2 -C 1 

and that (A^) ~ -q 2 /{a 2 ). Neglecting -In2 + 7E 

we easily find (Aq ^ ~ q 2 hiq. This result is only valid 

for q 2 /\ 2 <C 1 which corresponds, once plotted in Fig. 
^ to the solution for j — only. The case j > is 
solved using the fact that — (Inq)^ 1 <C 1: then the eigen- 
value equation is solved only when aj ~ vi.j which means 



lm(X") 

1 j = o 




FIG. 1: The first five (imaginary) eigenvalues Ay in units 
of i as a function of q for n = 0. 



(Af ) 2 ~-g 2 Mj forj>o. 

For n ^ instead, when q <C 1 we have 
qKi(q)/ Ko(q) ~ 2n. Then, since 2nJ n (x) ~ xJ n -i(x) + 
xJ n +i(x), we find J n _i(aj) = and, therefore, 

~ — Q 2 /v n -i,j- These asymptotic limits are com- 
pared with the actual solutions of the eigenvalue equation 
in Fig. □ Fig. Hand Fig. 

Note that in the region q <C 1, if the alternative nor- 
malization using Ap is selected, p shows only a weak log- 
arithmic dependence on the transverse beam size q in 
the case n = 0, j = Q and no dependence on q in the 
other cases. Looking at the slopes in the figures we can 
conclude that, with the use of A p in place of A p , p is, 
for realistic choices of /, much smaller than unity. The 
same applies when q — > 00: in this case A p ~ A p and p 
will also be small with respect to unity. As a result p, 
defined using A p , can be considered much smaller than 
unity in a wide range of parameters which justifies, at 
least in this particular situation, the assumptions used 
in the derivation of Eq. (|21|) . 

We can now write our final solution in the following 
form: 

p {n)(i ~ _ f Ej u nj J n { aj f)e x ^' z f < 1 

2 ( ' } [Z^^K^A 3 ' f>l ' 

and the coupling factors u n j are given by: 

K n {q)^d^J n {a 3 QU {n) {0 

J «( a j)dp" [aJ n +i(a)K n (q) - qK n+1 (q)J n (a)] p=x ( n) 

(91) 
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FIG. 2: The first five (imaginary) eigenvalues A^ in units 
of i as a function of q for n = 1. 



Ima (n) ) 




n = 2 



FIG. 3: The first five (imaginary) eigenvalues A^ in units 
of i as a function of q for n = 2. 



Moreover it must be finite at f = 0. We do not work out 
details, which can be found in 0], but we underline the 
fact that the anti-Laplace transform of E( n > calculated 
with this method coincides with our previous result Eq. 
(|90|l . Note that in the case the Green function is de- 
rived without the expansion in ty n j the final solution for 
E z is automatically valid, but still subjected to the as- 
sumptions on the validity of Jordan lemma; without the 
introduction of other assumptions on the system under 
study we can safely say that our result holds for the case 
of a cold beam only. 



C. Parabolic profile 

A parabolic transverse profile corresponds to the case 
S'o(r) = 1 — kff 2 . This is one of the few profiles for which 
the evolution problem can be solved analytically. The 
study of this situation offers, therefore, the possibility 
of crosschecking analytical and numerical results with or 
without the use of the semi-analytical method described 
in Section ITlTDl 

In the parabolic case g(f,p) inside the operator C is 
given by 



l(r,p) = 



„2 



1-fcf 



f < 1 
f > 1 



(92) 



Here we assume, strictly, F(P) — S(P). Solution for the 
homogeneous problem defined by C can be found in 0, 
since it is of relevance in FEL theory as well. We can use 
that solution in order to solve our eigenvalue problem, 
and to write the expressions for the eigenfunctions ^> n j 
to be inserted in Eq. I(63|) . Let us introduce the following 



notations: /i 2 — iDq 2 - 
e = {n + l)/2-n 2 /{4S). 



A 



5 2 



A 



(n) 



- ami, d 

After some calculation we find: 



_ _ J f n e~ Sf ' 2 / 2 iF x (e,n+ l,Sf 2 ) r<l 



f > 1 



where \F\ is the confluent hypergeometric function, and 
the eigenvalue equation analogous of Eq. 1)89(1 is now 



where a = ^/g(r,p)\ f<1 . 

At this point we should show that the expansion Eq. 
(|51|l is correct, so that the method used up to now can 
be rightfully applied. Yet we cannot do this. We assume 
this fact, instead, and we prove that this assumption is 
right both using numerical techniques in Section llVI and 
with the help of an alternative analytical technique here: 
in fact, interestingly enough, one can solve Eq. I|39|l also 
by finding directly a Green function without any partic- 
ular expansion simply imposing that the Green function 
obeys CG^ — for all f except f = f', where it must be 
continuous and such that its derivative is discontinuous 
(the difference of the left and right limit must equal 1 jf) . 



5K n (d) [2e(n+ l)" 1 iF x (e + 1, n + 2, S) 
- 1 F 1 {e,n+l,5)]+dK n+1 {d) 1 F 1 {e, n + 1, 5) = (94) 

Once more we should show that the expansion Eq. 1(5 1|) 
is correct. We assume this fact, and we present a cross- 
check of our result Eq. ((93(1 using numerical techniques 
in Section llVl 

D. Multilayer method approach 

An arbitrary gradient axisymmetric profile can be ap- 
proximated by means of a given number of stepped pro- 
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files, or layers, superimposed one to the other. This 
means that results in Section [ill Bl can be used to con- 
struct an algorithm to deal with any profile (see or 
|14| for more details and a comparison with the same 
technique in FEL physics). 

The normalized radius of the beam boundary is simply 
unity; let us divide the region < f < 1 into K layers, 
assuming that the beam current is constant within each 
layer. Within each layer k, according to Eq. l(2"9")l. the 
solution for the eigenfunction is of the form 

$L fc) = c k J n (fi k f) + d k N n (fj, k r) , (95) 

where (k — 1)K < f < k/K, c k and d k are constants, J n 
and N n are the Bessel functions of first and second kind 
of order n, and 

(4 = -q 2 (l - ibS k _ 1J2 ) . (96) 

Here Sfc- 1/2 = 5 , (r fe „ 1 / 2 ) and r k -i/2 = (k-l/2)/K. To 
avoid singularity of the eigenfunction at f = we should 
let d\ — 0. Then, the continuity conditions for the eigen- 
functions and its derivative at the boundaries between 
the layers allow one to find all the other coefficients. The 
continuity conditions can be expressed in matrix form in 
the following way: 



The two relations above can be also written in matrix 
form as: 




(102) 



where the coefficient b can be expressed in terms of the 
coefficient c\ by multiple use of Eq. I|97|l . Since c\ can 
be chosen arbitrarily, we may set c\ = 1 to obtain the 
following matrix equation: 



T K x T x _! x...xT 1 (^J=T(M=&IJJ . (103) 

The matrix T depends on the unknown quantity A. The 
other unknown quantity in Eq. \ 103(1 . the coefficient 6, 
can be easily excluded, thus giving the eigenvalue equa- 
tion: 

(T) n = (T) 21 , (104) 

which allows one to find the eigenvalue A. Then using 
Eq. (jHIl an d Eq. (|102(l it is possible to calculate the 
eigenfunction. 



/ c k+ i \ ( c k \ IV. APPLICATIONS AND 

( 4+i J =Tk [ d k J ' fc = 2 ' " ,K ~ 1 ' ( 97 ) EXEMPLIFICATIONS 

where the coefficients T k are given by (f k = k/K): a. Algorithm for numerical solution 



(Tfe)n = {ir/2)r k [[i k J n+ i(p k r k )N n (n k+ ir k ) 

-fJ-k+lJn(pkf k )N n+ i(fI k+ if k )] 

(T k )i2 = {ir/2)r k [ii k N n+1 (p, k f k )N n {p, k+1 f k ) 
-li k+1 N n {p, k f k )N n+1 (p, k+1 f k )] 

(T k ) 2 i = -(7r/2)f fc [n k J n+ i(n k r k )J n (p k+ ir k ) 
—^ k+ iJ n (^ k r k )J n+ i(fi k+ ir k )} 

{T k ) 2 2 = -(7r/2)ffe [n k N n+1 (n k r k )J n (p k+ ir k ) 
-H k+ iN n (n k r k )J n+1 (p k+ ir k )} 



(98) 



Eq. H29[) a l so dictates the form of the solution for the 
eigenfunction outside the beam r > 1, satisfying the con- 
dition of quadratic integrability: 



$„(f) = bK n (qf), Re(q) > 



(99) 



Then, continuity at the boundary, i.e. at f = 1 gives the 
following relations: 



c k J n (p-k) + d k N n (fj, k ) = bK n (q) (100) 



and 



H k c k J n+1 (p k ) + n k d k N n+1 (p k ) = bK n+1 (q) . (101) 



The results in Section ITTll constitute one of the few ex- 
isting solutions for the evolution problem of a system of 
particles and field. Yet, to derive it, we had to rely on 
several assumptions, among which that of a small pertur- 
bation, in order to get linearized Maxwell- Vlasov equa- 
tions. This is not too restrictive, since in practice one 
has often to deal with space-charge waves in the linear 
regime, but it would be interesting to provide a solution 
for the full problem. From this viewpoint, the only way 
to proceed is the development of some numerical code 
based on macroparticle approach capable to deal with 
the most generic problem. As a first, initial step towards 
this more ambitious goal we present here a numerical 
solution of the evolution equation in the case of an axis- 
symmetric beam, that we will cross-check with our main 
result, Eq. (|63|l . In order to build a numerical solution 
one may, in principle, use Eq. I(29|l . but it turns out more 
convenient to make use of Eq. 1(3 1|) . 

Eq. 1(31(1 can be specialized to the axis-symmetric case 
by integration of Eq. ((30(1 over the azimuthal angle 4> 
which leads to 



-iq 



dr 



£f(»)g— in<t> 



(105) 
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FIG. 4: E — Re(E z ) as a function of z and r. Here 9=1, 
7i = and the first five eigenfunctions have been used. A 
stepped transverse profile has been used. 





FIG. 6: E — Re(Ez) as a function of z and f. Here q = 
1, n = and the first five eigenfunctions have been used. 
A gaussian transverse profile proportional to e r A " ' with 
cr = 2.0 has been used. 
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FIG. 5: E — Re(E z ) as a function of z and r. Here q = 
1, n = and the first five eigenfunctions have been used. 
A parabolic transverse profile proportional to 1 — k\f 2 with 
k\ — 1.0 has been used. 



where 

(r ' r) ~\/„(<zf')^„(gr) f>? ' (iUbj 

I n being the modified Bessel functions of the first kind of 
order n. The equation for the n-th azimuthal harmonic 



FIG. 7: Comparison between analytical results (solid line) 
and numerical methods (circles). E — Re(E z ) is plotted as 
a function of f for several values of z. Here q = 1, n = 
and the first five eigenfunctions have been used. A stepped 
transverse profile has been used. 



of the field can be written from Eq. (|105fl as 

= ^ C df'f'j[ n) G^ . (107) 
Jo 

Therefore, under the assumption of a cold beam, Eq. i|31|) 
can be rewritten as 
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FIG. 8: Comparison between analytical results (solid line), 
multilayer approximation with 15 layers (dotted line) and nu- 
merical methods (circles). E = Re(E z ) is plotted as a func- 
tion of r for several values of z. Here q — 1, n = and the 
first five eigenfunctions have been used. A parabolic trans- 
verse profile proportional to 1 — k\f 2 with k\ — 1.0 has been 
used. 



Analytical 
Multilayer 




FIG. 9: Comparison between multilayer approximation with 
15 layers (solid line) and numerical methods (circles). E = 
Re(Ez) is plotted as a function of f for several values of z. 
Here q = 1, n = and the first five eigenfunctions have been 
used. A gaussian transverse profile proportional to e~ r ^ 2a -* 
with <r = 2.0 has been used. 



3i 



(n) _-("), „■"(") 



'id 



-q S / dz 
Jo 



(z'-z) / df'f'j^G^ 



(108) 



which can be easily transformed, by double differentia- 
tion with respect to z into the integro-differential equa- 
tion: 



d 2 j[ n) 
dz 2 



= -?S I dr'r'G^ 



(109) 



Eq. (|100|) is of course to be considered together with 
proper initial conditions for ji and its z-derivative at 
z = 0. The interval (0, 1) can be then divided into an 
arbitrary number of parts so that Eq. (|109fl is trans- 
formed in a system of the same number of 2nd order 
coupled differential equations. The possibility of trans- 
forming Eq. (|109f) in a system of 2nd order coupled 
differential equations explains the choice of starting, in 
this case, with Eq. l|31fl instead of Eq. I|29|l : in this way 
our system can be solved straightforwardly by means 
of numerical techniques. To do so we used a 4i/i-order 
Runge-Kutta integration method, which gave us the so- 
lution of the evolution problem in terms of the beam 
current. Then, using Eq. (|107fl we could get back E z 
and we compared obtained results with Eq. 163fl for dif- 
ferent choices of transverse profiles. The real field E z 
should be recovered, for any particular situation, passing 
to the dimensional quantity E z and, then, remembering 
E z = E z e % ^' + E* z e~ % ^: yet, all relevant information is 
included in Re(E z ). To give first a general idea of the 
obtained result we present, in Figs. 0][5]and|I)| Re(E z ) 
as a function of z and f in the case of stepped, gaussian 
and parabolic profile respectively, with parameters choice 
specified in the figure captions. In all cases the initial 
conditions are proportional to the transverse distribution 
functions (stepped, gaussian, parabolic), F(P) = 5(P) 
and n = 0. We consider only initial density modulation 
(i.e. di e = 0). Note that, in order to be consistent with 
the perturbation theory approach we should really choose 
did *C 1, since it is normalized to the bunch current den- 
sity. However using, for example, did = P with p« 1 
will simply multiply our results by an inessential factor 
p so, for simplicity, we chose p = 1. Note the oscilla- 
tory behavior in the z direction. Comparison with the 
Runge-Kutta integration program are shown in Figs. \7\ 
|S]and|5| In the parabolic case, both pure analytical solu- 
tion and solution with multilayer approximation method 
are present, while in the gaussian case only a solution 
with the multilayer method is possible, to be compared 
with the numerical Runge-Kutta result. This comparison 
shows that the assumption of the validity of Eq. i|51fl is 
correct in the parabolic case, and validates it once more 
for the stepped profile situation. 
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FIG. 10: The final result and the first five modes aid = 
e -r 2 /(2<r) with a = o.l at S = and for q = 1. The final 
result is the sum of the first fifty modes, here. 



FIG. 12: The final result and the first five modes aid = 
e -r 2 /V°) w ith cr = 0.1 at i = and for q = 100. The final 
result is the sum of the first fifty modes. Here n — 0. 
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FIG. 11: The final result and the first five modes aid = 
e-f- 2 /^) w ith cr = 0.1 at z = 10 and for q = 1. The final 
result is the sum of the first fifty modes, here. 



B. The role of the initial condition 

Here we present some further exemplification of the 
obtained results. 

In particular we are interested in investigating, in a few 
cases, what is the role of the initial condition in the final 
results, in order to develop some common sense regarding 
our analytical formulas. 



The main parameter in our system is the transverse 
beam extent q. However, intuitively, one can set two lim- 
iting initial conditions: one in which only a small part of 
the transverse section of the beam is modulated and the 
other in which all of it is modulated. Depending on the 
profile of the initial modulation, one can have excitation 
of many modes or only a few, while the conditions for ex- 
citation of a single mode have been discussed in Section 
IIII The transverse parameter q will fix the eigenvalue 
problem and, therefore, the oscillation wavelength (in the 
z direction) of the various modes. If q is smaller than or 
comparable to unity we expect to have appreciable dif- 
ferences in the eigenvalues, which lead to a quick (in z) 
change of the relative phases between different modes. As 
a result the initial shape of the fields will change pretty 
soon. On the other hand, when q is larger than unity, 
we will have all the eigenvalues converging to unity as in 
the one-dimensional case, which means that the relative 
phases between different modes will stay fixed for a much 
longer interval in z and the initial shape of the fields will 
not change during the evolution. 

To exemplify this situation we set up several calcu- 
lations using our analytical solutions to the initial value 
problem. In particular we considered two cases q — 1 and 
q = 100, and a radial stepped profile. Again, we consider 
only initial density modulation (i.e. a le = 0) and we 



study two subcases: in the first we set a\ c 



= p -t 2 /{2a) 



with a = 0.1 and in the second we put aid = 1. As said 
before, in order to be consistent with the perturbation 
theory approach we should really choose a±d *C 1, since 
it is normalized to the bunch current density but using, 
for example, aid — P with p <C 1 will simply multiply 
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FIG. 13: The final result and the first five modes Sid — 
e - f2 /( 2tr ) with a = 0.1 at z = 10 and for q = 100. The final 
result is the sum of the first fifty modes. Here n = 0. 



FIG. 14: The final result and the first five modes a\d = 1 at 
z = and for 5 = 1. The final result is the sum of the first 
fifty modes. Here n = 0. 



our results by an inessential factor p so, for simplicity, 
we chose, again, p = 1. Moreover we considered the az- 
imuthal harmonic n — 0. 

Figures ITUI and ITT1 present the first five eigenfunctions 
(with relative weights and phases) and the sum of the 
first fifty (i.e. the final result) for aid — e ~ r ^ 2<T ' ) with 
(7 = 0.1 at z = and at z — 10 and for 0=1. As one 
can see, the relative phases have changed and the shape 
of the total field, our final result, has also changed with 
z. 

On the contrary, Figs. IT21 and present the first five 
eigenfunctions (with relative weights and phases) and the 
sum of the first fifty (i.e. the final result) for Sid = 
e -f 2 /(2a) with cr ^ o.l at z = and at z = 10, and 
for q = 100. Here the relative phases have almost not 
changed and the shape of the total field, our final result, 
is also remained unvaried (of course one must account for 
the fact that the system is undergoing plasma oscillation, 
so the shape, and not the field magnitude, is what is 
important here). 

For comparison, it is interesting to plot analogous fig- 
ures for the second situation, that is Sid = 1- Figs. ITU 
and 1 151 depict the situation for q = 1 at z = and z = 10 
respectively. The way the phases behave is similar to 
what has been seen before, i.e. there is a rapid change 
in the relative phases between the modes, but now it is 
more difficult to see from the plots because, in contrast 
with the case of aid = e^ 7 "^ 2 ^ already the first mode 
is sufficient to fit the initial conditions relatively well so 
that the field shape is almost unchanged. In Figs. 1 161 and 
1171 we plot, instead, the case q — 100 always as z = and 
z = 10 respectively, and with aid = 1- Here it is easy to 
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FIG. 15: The final result and the first five modes a\d — 1 at 
z = 10 and for q = 100. The final result is the sum of the 
first fifty modes. Here n — 0. 



see, once more, that the relative phases between different 
modes are almost unchanged. What is of interest in this 
latter set of four pictures is the way different modes are 
working together to satisfy initial conditions, in compari- 
son with the way they mix in the case for aid — e~ r ^ 2<T - 1 : 
note, in particular, how the first mode is almost enough 
to satisfy the initial conditions (Figs. 1141 and lTSIl . while 
in Figs. El and ^2 it is almost completely suppressed. 
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FIG. 16: The final result and the first five modes aid = 1 at 
2 = and for q = 100. The final result is the sum of the first 
fifty modes. Here n = 0. 



FIG. 18: The final result and the first five modes when the 
initial condition is set in order to excite the third mode alone. 
The final result is the sum of the first fifty modes. Here n = 0. 
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FIG. 17: The final result and the first five modes aid = 1 at 
z — 10 and for q — 100. The final result is the sum of the 
first fifty modes. Here n = 0. 



FIG. 19: The final result and the first five modes when the 
initial condition is set in order to excite the third mode alone. 
The final result is the sum of the first fifty modes. Here n = 0. 



This should not be a surprise, considering that we may 
actually select a single mode by fixing appropriate initial 
conditions as described in Eq. J7TJ and Eq. J72J- For 
instance, if we fix d le = and we want to excite only the 
jth mode for a given value of the azimuthal harmonic n 
in the case So = 1, then, according to Eq. i|72|) and Eq. 



(19011 we must set (modulus a constant factor): 

(110) 

where ctj is defined in Eq. I|89fl • For example if we choose 
n = and j = 2 (third mode) for q = 1 we obtain the 
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results presented in Fig. ^]at z = and in Fig. 1191 at 

z = 10. As it can be seen by inspection only the third 
mode is excited and evolves. Our condition Eq. H110JI set 
strictly to zero the contributions of all the modes with 
j =/= 2. Of course, in practice, actual data plotted in 
the figures show finite contributions of the other modes 
ascribed to the finite accuracy of our computations: to 
be precise, the difference between the final result (sum of 
the first fifty modes) and the third mode alone was found 
to be on the fourth significative digit. 

V. CONCLUSIONS 

In this paper paper we presented one of the few self- 
consistent analytical solutions for a system of charged 
particles under the action of their own electromagnetic 
fields. Namely, we considered a relativistic electron beam 
under the action of space-charge at given initial con- 
ditions for energy and density modulation and we de- 
veloped a fully analytical, three-dimensional theory of 
plasma oscillations in the direction of the beam motion. 

We used the assumption of a small modulation so that 
we could investigate the system behavior in terms of a 
linearized Vlasov equation coupled with Maxwell equa- 
tion, under the assumption that field retardation effects 
can be neglected. Then we introduced normalized quan- 
tities according to similarity techniques and we provided 
two equivalent presentations for the evolution problem 
in terms of a integrodifferential equation for the electric 
field and of a integral equation for the beam current. 

The integrodifferential equation for the fields was par- 
ticularly suited to be solved with the help of Laplace 
transform techniques: we did so in all generality and 
we discussed the mathematical difficulties involved in 
the general treatment, namely the assumption of a well- 
behaved differential operator allowing eigenmodes expan- 
sion of the Green function and the problem of the ana- 
lytic continuation of the Laplace transform of the field to 
all the complex plane (isolated singularities excluded), in 
relation with the application of the Fourier-Mellin inte- 
gral to antitransform E. Our considerations led us to 
restrict our attention to the cold beam case. We special- 
ized the general method to the important cases of stepped 
and parabolic transverse profiles, which are among the 
few analytically solvable situations. In particular, the 



stepped profile case could be used to develop a semi- 
analytical technique to solve the evolution problem for 
the field using an arbitrary transverse shape. 

We tested our results by discussing the limit for the 
1-D theory (q — » oo). We also developed an algorithm 
able to solve the evolution problem in terms of the beam 
currents. The integral equation for the currents could be 
easily approximated to a system of second order ordinary 
differential equations which could be solved by means of 
numerical Runge-Kutta integration method. Once the 
solution for the current was known we recovered the elec- 
tric field evolution by integration of the current with a 
suitable Green function. Numerical and analytical or 
semi- analytical solutions for the fields were then com- 
pared and gave a perfect agreement. In this way we 
could state that the assumption of the correctness of the 
eigenmodes expansion for the Green function has been 
proved, for some particular profiles, by means of numer- 
ical crosschecks (in the stepped profile case, by means of 
alternative analytical techniques too). 

Finally we exemplified the role of the initial condition, 
which we have seen to control the way one ore more 
modes interact together to give the final result. In partic- 
ular we have shown how to build up initial conditions in 
such a way that a single mode is excited and propagates 
through. We checked our prescription by setting up par- 
ticular initial conditions and looking at the propagation 
of various eigenmodes. 

In conclusion we proposed, checked and analyzed, both 
from physical and mathematical viewpoint, a theory of 
space-charge waves on gradient-profile relativistic elec- 
tron beams. This work is of fundamental importance, 
since it is one of the few known analytical solution to 
evolution problems for systems of particles and fields. In 
particular, today, it is of great relevance in the physics of 
FEL and high-brightness linear particle accelerators. 
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